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INTRODUCTION 


Boston  College,  in  its  investigation  of  the  seismic 
effects  of  rocket  launchings  on  structures  in  the  immediate 
environment,  is  analyzing  data  taken  during  a  Titan  IIT-D 
launch  at  Vandenberg  AFB  in  March  1979-  The  work  reported 
herein  deals  with  the  spectral  form  of  the  induced  surface 
pressure.  A  theoretical  model  is  fitted  to  the  power  spectral 
density  of  the  observed  pressure  at  various  times  during 
launch.  The  results  are  presented  along  with  tests  of  their 
validity. 
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CHARACTERIZING  PRESSURE  SPECTRA 


'T’he  data  are  contained  on  a  file  (LAUNC2.MNR)  which 
begins  1.64  seconds  before  ignition  and  continues  for 
61.44  seconds  at  100  samples  per  second  on  each  of  16 
channels.  The  array  of  instruments  corresponding  to  these 
channels  is  shown  in  Fig.  1.  For  this  report  only  the 
pressure  sensor  on  channel  10  which  is  050  feet  from  the 
launch  point  was  used. 

To  obtain  estimates  of  the  power  spectral  density  at 

a  given  time  after  ignition  256  samples  were  read  in 

beginning  at  that  time.  These  samples  had  been  digitized 

at  204.8  counts  per  volt.  The  power  spectra  were  computed 

by  the  periodogram  technique  employing  on  FFT  algorithm. 

1  N  z  N 

This  periodogram  is  such  that  ^  E  x.=  T.  S,  . 

Ni=l  1  k=l  K 

It  was  found  that  the  noise  caused  by  the  system  and 
ambient  pressure  fluctuations  was  white  and  assumed  to  be 
additive  above  the  quantization  level.  It  was  estimated 
by  averaging  20  periodograms  starting  1  minute  before  launch 

_C  •p 

and  found  to  average  9.37(10  )  volts'  per  cell  (Fig.  2). 

This  was  subtracted  from  the  rocket  spectra  to  yield  our 
best  estimate  of  rocket  induced  surface  pressure  observed 
through  our  system.  To  convert  these  spectra  to  the  spectra 
of  the  pressure  appearing  at  the  input,  they  were  divided 
by  the  squared  magnitude  of  the  system  response  (Fig.  3). 
That  is  =  k  .  These  were  scaled  by  1/af  at  Ohz  and 


the  Nyquist  frequency  and  2/af  elsewhere  so  the  resultant 


2 


power  spectral  density  when  integrated  by  trapezoidal  rule 
from  Ohz  to  the  Nyquist  frequency  equals  the  mean  square  of 
the  input.  However  since  the  value  at  Ohz  is  dependent  on 
amplifier  drift  its  value  is  discarded  and  to  reduce  contamina¬ 
tion  from  errors  in  our  estimate  of  system  behavior  near  the 
Nyquist  frequency  only  frequencies  less  than  40hz  were 
considered. 

From  physical  considerations  (1)  and  experimental  studies 
‘2)  it  is  believed  that  the  power  spectral  density  of  the 
surface  pressure  caused  by  undeflected  chemical  rocket  plumes 
is  of  the  form: 


,  P  „  W  -2 
P(w)=l  -***  +-0} 

V  ‘  TT  WQ  XW0  W  ' 


4  P0  f  f0 
or  P(f)=i  ^ 


-2 


It  was  desired  to  obtain  values  of  PQ  and  fQ  which  minimized 

the  sum  of  the  squared  differences  between  the  observed  power 

spectral  density  and  the  theoretical  as  described  by  the  above 

equation.  To  find  the  sum  of  squared  errors  over  the  range 

of  the  power  spectral  density: 

129 

E=  I  (P'-P(fk))2 
k=l  k 

where  P£  is  our  estimate  of  the  k*"*1  value  of  the  power  spectral 
density  and  fk=100 (k-1) /256  which  is  the  frequency  represented 
by  the  kth  values  of  the  power  spectral  density. 

For  reasons  already  stated  the  full  range  of  the  power 
spectral  density  was  not  to  be  used  so  the  limits  of  the 
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summation  were  changed  accordingly. 


103  - 

E=  (P'-P(f.)p 
k=2  k 


Rather  than  finding  both  Pg  and  fg  by  trial  and  error 
the  derivative  of  E  with  respect  to  Pg  is  taken  and  set  to 

zero  that  is  where  the  extremal  would  be  found. 

103  4Pn  f.  f  -2  .  f.  fn-2 

||=0=  Z  2  (P ' 0}  )  (-If  {  *+  0}  ) 

dP0  k=2  k  f0  0  fk  0  f0  f k 


irf 

P  =~ 4 
0 


103  £  f„  -2 

£  p  •  (-  lr_8 — 

o  k:2  k  Vfk 

103  f.  fn  7 

£ 

k=2  0  "k 


So  for  any  value  of  fQ  the  value  of  PR  which  insures  minimum 
error  in  the  least  square  sense  is  uniquely  determined. 

To  find  the  value  of  fQ  which  yields  the  smallest  E 
(the  'best'  fg)  an  iteration  scheme  was  devised.  The  initial 
search  interval  has  f£  and  as  its  endpoints  in  the  belief 

that  fg  lies  between  them.  Five  trial  values  of  fg  are  con¬ 
sidered  starting  with  the  low  endpoint  of  the  search  interval 
and  increasing  in  equal  steps  to  the  high  endpoint.  For  each 
trial  value  of  fg,  Pg  and  E  are  computed  by  the  formulae  above 
and  the  best  fg  selected.  A  new  search  interval  is  defined 
with  endpoints  equal  to  the  trial  values  adjacent  to  the  best 
fg  during  previous  search  and  a  new  estimate  of  the  best  fg 
found.  The  process  is  repeated  until  fQ  is  determined  to 
within  the  limits  of  the  computer's  accuracy.  In  this  case 
single  precision  yields  about  seven  significant  digits. 


The  resultant  fg  and  Pg  were  used  to  qenerate  plots 
which  show  power  spectral  density  based  on  observed  values 
and  theoretical  values  vs.  normalized  frequency,  f/fg. 

Table  I  shows  fg  and  PQ  for  each  segment. 

In  order  to  compare  average  observed  power  with  theoretical 
power  the  theoretical  power  spectral  density  was  integrated 
from  to  00 . 

pt=V  r-J  (f  *r-]  2,3£ 

C  .1  Xg  Xg  X 

A  c  r°°  f2 

=?  Vo  '  7J— 2.df 

(f  +fg) 

_4P0f0  1— S5 — I — 5*  T  +?5  -f  — 5 — J  df  (3] 

*  1  <f +f0}-i  -  (f  +f0} 


2P  f  2P  f 

. 0E0  1  ZtVo,  n 

- - -  1  — - -  df—  — — — ( — 2 — )  — 2P 

"  (f 2+f q)  "  £0  0 


The  observed  power  (Ax  E  P')  and  theoretical  Dower  are  contained 

k=l  K 

in  Table  I. 


VALIDATION 

Although  the  theoretical  curves  do  not  appear  to  closely 
fit  the  observed  spectra  it  is  known  that  the  use  of  the 
periodogram  with  a  large  number  of  samples  to  estimate  the 
spectrum  produces  results  which  fluctuate  wildly  (4).  It 
was  desired  to  test  whether  these  fluctuations  fall  within 


expected  limits  when  the  theoretical  power  spectral  density 
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is  assumed  true. 


Oppenheim  and  Schafer  (4)  show  the  development  of  an 

expression  for  the  variance  of  the  estimated  's  for  a 

white  Gaussian  process  and  Hinich  and  Clay  (5)  state  that 

the  result  is  a  good  approximation  for  a  wide  variety  of 

2 

random  processes.  They  state  that  the  variance  (o  )  of  the 
spectrum  is  approximately  equal  to  its  magnitude  squared. 


v  (p ' -p * } 2=p 1 ^ 
°k  Mm-1  Pk  k.  k 


or  %  _  i  "  (p; <m)-p/)2 

2  m  _ _ _ _ _ 

m=l  pi  2 

1  103  ak 

A  figure  of  merit  is  defined  as  Z  - =- 

k=2 


which  'with  M=1  and  assumed  equal  to  P(£k)  becomes 

i  x?3  (p:-p(fk))2 
102  k=2  --2~- - 

p2(f  ) 


This  quantity  which  should  be  approximately  1  is  contained  in 


Table  I.  In  addition  as  has  been  stated  in  (5),  for  large  N 
P 1 

2 _ k_  follows  a  chi-squared  distribution  with  2  degrees  of 

P(fk) 

freedom. 

As  an  additional  test  of  the  validity  of  the  fit  this 

9T>  1 

criteria  was  used.  For  each  value  of  A,  =*  *  k  its  cumulative 

ft^> 

relative  frequency  was  computed. 

a^=  estimated  probability  that  Ak<x. 
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acceptable  fits  simulated  data  was  used.  A  random  number 


subroutine  was  used  to  produce  normally  distributed  variables. 
For  each  group  of  256  numbers  the  FFT  was  taken  and  the 
magnitude  multiplied  by  the  square  root  of  P(f)  for  a  set 
value  of  Pq  and  fg.  "’he  inverse  transform  was  taken  and 
the  result  multiplied  by  a  decaying  exponential  which  modelled 
the  envelope  of  the  rocket  data.  Two  hundred  similarly 
produced  groups  of  256  variables  were  fitted  with  the 
theoretical  power  spectral  density  curve  and  statistics  of 
the  figures  of  merit  accumulated. 

The  figures  of  merit  for  this  simulation  were  found  to 
have  a  mean  of  1.07  and  a  standard  deviation  of  .35.  The 
minimum  value  obtained  was  .57  and  the  maximum  was  2.98. 

It  was  decided  to  accept  the  fitted  data  if  its  figure  of 
merit  fell  within  these  extrema. 

The  first  fit,  which  began  3.83  seconds  after  ignition, 
produced  a  figure  of  merit  of  21*1.6.  This  fit  is  rejected. 

Tt  is  believed  that  this  early  in  the  launch  the  plume  was 
not  undeflected  and  the  theoretical  curve  does  not  apply. 
Figures  of  merit  for  subsequent  fits  fell  within  extrema 
criterion  and  these  fits  are  accepted. 

It  was  found  that  when  f  was  low  (eaual  to  2.8  hz) 

o 

the  chi-squared  test  plot,  while  a  straight  line,  lay  above 
the  A. =u  line,  as  happened  on  a  number  of  the  fitted  segments, 
and  the  estimated  value  of  f  was  high  (3.2hz).  Evidently 
out  estimates  while  acceptable  according  to  the  figure  of 
merit  criterion  are  biased  toward  the  high  frequencies. 
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APPENDIX 

Listincr  of  computer  program  follows 


FORTRAN  IV 


VU? 


Fri  ,J  9  "Nov -7  9  00:00: '.)!) 


A  i. 


c 


c 

c 

c 

c 


OV-  3 

0,104 

.5  0.15 
*  1 0,*.  f, 
3  or-  7 


ODD  9 
0010 
noil 
no  12 
no  n 
0314 
0015 
no  15 


c 

is 


]  5 


0017 
00 10 

r 

0019 
00  20 

C 

0071 
0  0  2  2 
00  22 

t'O  24  70 

C 

n,  2  5 
0025 

?■  I?" 

0029 

C 

0  0  29 


PROGR.V*  POWER 4 

THIS  IF  POWER 2  WITH  MUCH  OF  THE  OUTPUT  LEFT  OUT. 

LINK  THIf.  WITH  PLLIB 3  FOR  PLOTTING. 

LINK  WITH  MIRC  FOR  HISTOGRAM  AND  ASSOCIATED  DATA  TABLES, 
DATMIN  AND  DATMAX  FUNCTIONS,  AND  SEARCH  SUBROUTINE. 
LINK  WITH  FFTSUB  FOR  FAST  FOURIER  TRANSFORM  SUBROUTINES. 


DIMENSION  WORK  (102),  PTHEOR  (900  )  ,  ASCALF.  (9),  IN  (20),  DUMMY  (255'  , 
CONV  (109 )  ,  NSKIP  (12),  RESID(102),  AJUNK(75),  PMlJLT(/,> 
STATS  (102),  THSTAT  (102)  ,  DAT  (  5 )  ,  VARPLT  ( l  ('•!  '  ,  SSL  PLT  ( > 
LOGICAL*!  LABEL  (5) 


COMMON  PI,  W ( 3  0  0 ) ,  PPRIME  (255 ) 

DATA  PI/9 . 14 15925535999/ 

DATA  ASCALF./ 175. 9,131.4,204.8,171.4,209.9,153.9,  152.9,204 
DATA  NSKIP/0 , 235 , 39 1 , 54  7 , 7g  3 , 859 , 1 0 1  5 , 1  1  7  1  ,  L 1 27  ,  1 4  8  3  ,  1  5  ’0 
DATA  P MULT/ .  57 5 , 1 . ,  2 . ,  2 . 77  3 / 


20  3 . 2/ 

7  9  r.  , 


TYPE  10 

FORMAT  (  '  TYPE  DATA  FILE  NAME:’  /) 

CALL  ASSIGN  (1,  '  -1,  ' RDO ’ ) 

DEFINE  FILE  1  (5144,  20,  U,  INDEX) 

CALL  ASSIGN (2 , ' DK : VF2PRS . DAT ’ ,13, 'ROO' ) 

DEFINE  FILE  2  (2 57 , 6, U , JNDEX ) 

TYPE  15 

FORMAT  ('  TYPE  FILE  NAME  FOR  OUTPUT  OF  RESIDUALS  AND  THEORETICAL’ 
+  /  •  VALUES:'  /) 

CALL  ASSIGN (2, '  ' ,-l ) 

DEFINE  FILE  2  (102 , 50 , U , KNDEX) 

DELTAF  =  100. /? 55. 

RNOYZ  =  . 935E-7 


.INDEX  =  l 
DO  20  1=1,103 

READ  (2  'JNDEX ) XXX , CONV ( I ) , YYY 
JNDEX  =  JNDEX  +  1 

TYPE*, 'CHANNEL  NUMBER  =  ?' 
•ACCEPT*,  1  CHAN 
TYPE *,' I  START  =  ?' 

ACCEPT*, ISTART 


CALL  DATE (DAT) 


00  30 

n 

DO  1003  ISKIP  =  ISTART, 12 

00  31 

PRINT  25,  DAT 

00  2? 

25 

FORMAT  ('  PROGRAM  POWER 4  '  1 0X , 3A 4  ,  // 

+  '  FILES  USED: ' l 3X 'LAUNC2.MNR  INPUT'/ 

+  25X 'VF2PRS.DAT  INPUT’  /  2 5X 'DATA. LC ? 

OUTPUT'  // 

0  <1  3  2 

PRINT  30,  ICHAN,  NSKIP  ( I SK  IP )  ,  RNOYZ /DELTAF 

00  2.1 

30 

FORMAT  ( ' N *****CHANNEL ' ,  T3,  10X,  T5,  '  OBSERVATIONS 

SKIPPED.’ 

33 


u  c..'  c.'  (.■  v  \r  v.-  c.j  r-  ; .  <  •  u  c> 


FORTRAN  [V 


V02.04 


Fri  09-Nov-79  00:00:00 


PAGE  002 


0.)  5  5 
vin3ri 
Of-P” 
o,;38 

0o 39  50 


o'J/JO 


r 

r 


q 


C 


■')S; 
OO.'P 
'JO A  ' 


0  ,! 


'»  * 


r 


4 


NOISE  SUBTRACTED  =',1PE11 


A) 


READ  IN  NEW  OBSERVED  DATA 


INDEX  =  N5K1P  (ISKIP)  +  1 
DO  50  I -1 ,255 
DUMMY  (I)  =  0. 

READ  (l  ’INDEX) IN 

PPRTME(I)  =  I N  (I CHAN  ) /ASCALE  (I CHAN -7  ) 

SUBROUTINE  FFTFP  CONVERTS  ARRAY  DATA  INTO  SPECTRUM 
CALL  FFTFP  (PPRIME , DUMMY ,256,2,0) 

CONVERT  SPECTRUM  INTO  PSD  (PPRIME  ARRAY).  DROP  FIRST  POINT. 
GENERATE  FREQUENCY  AND  CONVERT  TO  OMEGA  (W  ARRAY). 

•"ME  END  OF  THE  LOOP  IS  102  INSTEAD  OF  128,  SO  THAT  THE  MAXIMUM 
FREQUENCY  IS  NOW  40  HZ,  INSTEAD  OF  50  HZ. 

P3UM  =  0. 

Dn  50  1=1,102 

PPRIME  (I  )=2  .*  (PPRIME  (l+l  )-RNOYZ)  /  (QELTAF*CQNV  (T+l  )*CONV  (T+l  ) ) 
IF(PPRTME(I)  .LE.  0.)PRINT  58  ,  FLOAT  (T  )  *DELTAF,  PPR IME  ( C  ) 

FORMAT ( '  *****  WARN I NG  AT',  F5.1,'  HZ  PSI**2/HZ  =',G11.4) 

PS U M  =  PSUM  +  PPRIME  (I) 

W(T)  =  FLOAT  (I)  *  DELTAF 

I F ( I  START  . ME.  ISKIP)  GOTO  75 
KNDEX  =  1 
DO  70  1=1,102 

READ  (3 ’KNDEX) AJUNK 
AvIUNK  (2  5)  =  W(I  ) 

KNDEX  =  KNDEX  -  I 
WRITE  (3 'KNDEX) AJUNK 

FTND  THE  TRIAL  FREQUENCY  THAT  GENERATES  THEORETICAL  VALUES  THAT 
BEST  MATCH  THE  OBSERVED  VALUES  (CRITERION:  LEAST  SUM  OF  SQUARED 
ERRORS).  STORE  ANSWER  IN  WO.  STORE  ASSOCIATED  BEST  P  IN  P0. 


''ALL 
n.i  = 


SEARCH (WO , SSE , A , 0 ,W ( I ) ,W ( 1 92 ) , 0. , 7 ,WORK , PPR IME ,102,0,0- 
7. , 5. , ’TRIAL  FREQUENCY' , 1 5 , VARPLT, SSEPLT, 100) 
PNOT (WO) 


USING  THE  ESTIMATES  OF  P  AND  W,  CALCULATE  THE  CURVE  {PTHEOR) , 
THEN  CALCULATE  THE  RESIDUALS  AND  THE  FIGURE  OF  MERIT  (DOLLY) 
CALL  FUNCr (WO, PTHEOR) 

RESBAR  =  0. 

DOLLY  =  0. 

DO  550  1=1,1(52 

RESID(I)  =  PPR  IMF,  (I )  -  PTHEOR  (T  ) 

HEGBAR  =  RESBAR  f  ABS  (RESTD  (I ) ) 

HOLLY  =  DOLLY  F  (RESID (I ) /PTHEOR  (I ) ) **2 
RESBAR  =  RF.SBAR/  IW2 . 
h  •:  i  Y  =  noi.LY/lf<2. 
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C 

r 

C 

30  58 

0359  575 


PRINTOUT  RESULTS 

PRINT  575,  WO,  PO,  SSE,  DELTAF*PSUM,  2.*P0,  RESBAR,  HOLLY 
FORMAT  (///'  F-NOT  ESTIMATE  (HZ )  *  ,  T32 ,  1PF.35.7  / 

*  •  P  ESTIMATE ' ,  T.32,  1PE15.7  / 

+  !  SUM  OF  SQUARED  ERRORS  (SSE)’,  T32,  3  PE  IS.  7  / 

+  '  TOTAL  OBSERVED  POWER',  T32,  3  PE  15. 7  / 

v-  *  total  theoretical  power*,  T32,  ipeis.7  / 

5-  '  AVERAGE  ABSOLUTE  RESIDUAL’,  T’7,  1PE35.7  / 

»-  ’  FIGURE  OF  MERIT'  ,T32 , 1PE 1  5 . 7  ,  ////  ) 


337s; 


0971 
007  2 
907  3 
007  4 
007  5 
9075 

0  077 
0  0  ?  8 
007  9 

0030 
O  OBI 


PRINTOUT  RESIDUALS  AND  THEORETICAL  VALUES  WITH  FREQUENCY. 
STORE  THEM  IN  DATA  FILE  3. 


C  THE  FOLLOWING  PRINTOUT  IS  OPTIONAL 

CC  PRINT  593 


590 


FORMAT ( ’  FREQ 
+  1  (HZ  1 


THEORETICAL',  2{  39 X  'FREQ 
PST  **2 /HZ '  oX  'RESIDUAL',  7(1  AX 


THEORETICAL 
'  (HE '  P 


*-  6X  ’RESIDUAL')) 

CC  DO  593  1=1,34 

CC60U  PRINT  613,  (  (W{J)  ,  PTHEOR  (.J  )  ,RESID  (J  ) )  ,.7  =  I  ,  (I +58)  ,  34) 

610  FORMAT  (IX, F5. 2,  2  (5X, 1PE9.2) ,0p ,  2{15X,  F5.2,  2(SX,  1PE9,?) 

PRINT  625 
625  FORM  AT  (M*) 

KNDEX  =  1 


650 


DO  550  I=i , 102 

READ  (3  'KNOEXJAJUNK 
LDEX  =  TSKIP*2  -  1 
AJUN'K  (LDEX)  =  PTHEOR  (I  ) 

A  JUNK  (LDEX  f  1)  =  RESID(I) 
KNDEX  =  KNDEX  -  1 
WRITE  (3 'KNDEX) AJUNK 


C  OPTIONAL  PLOT 

C 

Q 

C  PLOT  PPRIME  (OBSER'*ED)  VS  PTHEOR  (THEORETICAL) 

C  INCLUDE  A  COMPARISON  LINES: 

C  1)  {  .575  *  PTHEOR)  VS  PTHEOR 

C  2)  (1.000  *  PTHEOR)  VS  PTHEOR 

C  3)  (2.O00  *  PTHEOR)  VS  PTHEOR 

C  4)  (2.773  *  PTHEOR)  VS  PTHEOR 


l  J  }  H  / 

;HJS  3 
U  8  * 
.0035 
0  35  5 
!'r;  =  7 
v  U  88 
■1939 


r 


PLOTTING  INFO: 

PTMIN  =  AMIN1 (PTHEOR (13 , PTHEOR (1G2)) 
PPM  IN  =  DATMIN  (PPRIME, 102) 

YMIN  =  AMIN  I (PPMIN , , 5*PTM IN ) 

PTMAX  =  D ATM AX (PTHEOR, 132) 

PPM AX  =  OATMAX  (PPRIME,  3  02  ) 

YMAX  =  AMAX1  (PPMAX , 2. 77  3 *PTMAX ) 
DELTAY  =  (YMAX  -  YMIN}/6 . 

DELTA X  =  (PTiMAX  -  PTMIN )/7. 
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l’r  i  iio-Nnv-vn  on*  on-,  no 


pack  -o-i 


00"P 


ooro  vis  mm.ot  s s  skipped. 


OOP  i 
;PPi;i 
o.jt)  : 
(S’.  '1 
OOP' 

ooo-. 

Poo  / 
ooo  . 

PO‘1'1 
0  1  00 


o  i  o  i 
S'  I  O'1 


’  1  0 


-I  1  s 

!  lit 


PALL  PLOTS  (0  ) 

CUL  PLOT  ( 1  ,  ,  -  ri .  ,  J) 

CALL  AX  I  S  1  (0 ,  ,  (1,  ,  1  OBSERVED  PS.  i  * *'>  -01'.'  1  ,  l  H  ,  f*.  ,  "0  .  ,  YM  !  N  ,  Hi  LTAY  ) 

PALL  AXISl  0,  VIlKORKTl  CAL  PS  l  *f  7 /W.  '  , -7.1  ,  ' .  ,  0 .  ,  PT  *1 !  N  ,  OELTAX  ) 

PM,L  AX  I  S  1  (7.  ,0.  ,  '  OBSERVED  PS  1**2/117.  '  ,  -10,6.  , '»S  .  ,  YM  1  N  ,  01  (  TAY  ) 
PALL  LINE  (»TIIEOR#PTMlN,l'PLTAX,PP!UMEf  YM I N , DELRAY ,  U'7  ,  1  ,  -  1  ,  '  '*  *  ) 
PLOT  THE  I  COM  PAR  1  SON  LINKS 
iv ‘  !U)  1-1,4 

PAUL  PLOT  (0  .  ,  (PMUl.T  (I  )  *PTM  IN  -  YM  I  N  )  /DELTAY  ,  3  ) 

('ALL  PLOT  (7  .  ,  (PMULT  ( I  )  *PTMAX  -  YM  IN  )  /DELTAY  ,  ?, ) 

PRINT  A? S 

CALCULATE  AND  PLOT  TEST  STATISTICS  Vc>  FREQUENCY 

o.t  ;?  n  i  *  t ,  i  o 

STATS  ( I  )  •  7  .  *  PPR  1  ME  { V  )  /PTUF.OR  ( I  ) 


('  OPTIONAI  PLOT 

P 

'  i  s||1  P.OTO  7 !» S 

pi  !H  PALL  DOPLOT  (7 A. ,W#STATS, 107 FREQUENCY  (MZ)  ’  ,  M  , 

•  ’OBSERVED  TEST  STAT I  ST  IPS  '  ,  74  ,  -  l  ,  '  *  '  ,  I  ,  1  ) 

0  1  OR  PRINT  "7 S 

P 

r  .'PTiOvlAL  III  STOP.R  AM  OF  TEST  STATISTICS 

PP  PALI,  i!  I  ST  (7  ,  ,  A.  .  0 ,  '  OBSERVED  TEST  STAT  l  ST  ICS'  ,24,  STATS  ,  1  >7,0, 1,U) 

>*  lRDER  DATA  (NOT  NECESSARY  IF  HIST  WAS  USED): 


l|0  >  '  7S 

DO  7  IS  1  ,  111) 

7  10  ’ 

MIN  »  1 

0  |  OP 

DO  /  10  .1  -\  ,  \  07 

1  1  O'! 

1  F (STATS.  (.11  ,LT.  STATS 

0  1  11  7  10 

CONTINUE 

0  1  1  ’ 

S  -  STATS*  (1  ) 

0  !  !  1 

STATS.  (1)  *  STATS  (MIN) 

■HI  4  M>' 

\ 

STATS  (MIN)  •*  S 

PA  I  (’111  ATE  THEORETICAL  TEST  STAT I  STt\  S  BASED  ON  1NVERS>E  IF 
Ml  SOU A HE  DISTRIBUTION. 

THIS  IS  DONE  FOR  THE  FIRST  101  OF  Tllr,  107  POINTS;  THE 
Till- ORET  ORAL  VALUE  FOR  THE  LAST  POINT  IS  INFINITE. 


I'l'. 
•  !  I  > 
0  I  i  / 


•1  ■> 


IV)  /•  IP  I-!,IO| 

OB ’.PUB  *  l  107. 

nisTATU  i  -7 .  *  a i, or,  {, ,  onsPRB) 


I-l  'T  THEORETICAL  VS  OBSERVED  VALUES, 
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9116 
0119 
01  20 
0121 
0  1  22 
0123 
0  j  2  4 


01 2  6 
9127 
0123 


0129 
'U39 
0131 
013  2 
013  3 
0  l  2  1 
0135 
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0135 
0137 
0138 
0139 
0 1  4  0 
0  14  1 


0144 

0  l  4  5 


0146 
014  7 
0148 
0149 
0153 
0151 
0152 
0153 


015  5 
0155 
0157 
0158 


0,,10) 


DELTAY  =  AMAX1  (STATS  (1  02)  , THSTAT  (101  )  )/6. 

DELTA*  =  THSTAT  (101 )/7. 

CALL  PLOTS  '6) 

CALL  RED 

CALL  PLOT (1 . f -6. 5,-3) 

CALL  AXI.S1  (0.  ,0.  , ’OBSERVED  TEST  STATISTICS  '  ,  ?4 , 6.  ,  90  .  ,  U.  ,  DELTA  Y  ' 
CALL  AXIS1  (0.  ,0.  , ’THEORETICAL  TEST  STATISTICS ', -27 , 7  .,  9.  , 

DELTAX) 

CALL  AXIS1 (7. ,0. , ’OBSERVED  TEST  STATISTICS -24 , 6. , 90. , 0. , DELTAY ) 
CALL  LTNE (THSTAT, 0. , DELTAX, STATS , 0. , DELTAY , 101,1,-1,'*") 

CALL  PLOT  {  THSTAT  (1 )/DELTAX  ,  THSTAT (1 ) /DELTAY  ,3) 

CALL  PLOT  (  THSTAT  (101) /DELTAX  ,  THSTAT (101) /DELTAY  ,?) 

N5KTP  EQUALS  SECONDS*100  PASSED 
TIME  =  FLOAT  (NSKTP  {'(SKIP)  )/l00.  -  1.64 
ENCODE  (6 , 759,  LABEL  (1  )  )TIM£ 

FORMAT  { F  8 . 2 ) 

CALL  SYMBOL  (1 ., 6. , 0 , 'CHANNEL  10',0„,10) 

CALL  SYMBOL  {3. ,6. ,0, LABEL, 0. ,6) 

CALL  SYMBOL  (3 . 75 , 6. , 0 , ' S ECONDS  AFTER  IGNITION 0. , 22 ) 

PRINT  625 

THE  FOLLOWING  EXTENDS  THE  PTHEOR  AND  W  ARRAYS  SO  THAT  THE 
THEORETICAL  CURVE  IS  SYMMETRIC 

NPOINT  =  THE  NUMBER  OF  POINTS  NECESSARY  FOR  A  SYMMETRIC  CURVE. 

NPOINT  =  101 

INDEX1  =  102 

NPOINT  =  NPOINT  +  1 

W  ( N  POINT )  =  FLOAT (INDEX1 )*DELTAF 

PTHEOR  (NPOINT)  =  PCURVE  (W0 , P0 ,W (NPOINT ) ) 

INDEX!  =  TNDEX1  +  10 

TF(PTHEOR  (NPOINT)  .GT.  PTHEOR  (.1  ) )  GOTO  ^75 

NORMALIZE  W  BY  DIVIDING  BY  W0 
DO  300  T=l, NPOINT 
W ( I )  »  W ( I ) /W0 

LOG-LOG  PLOT  OF  BOTH  OBSERVED  AND  THEORETICAL  VALUES  VS 
FREQUENCY. 


PM  IN 
PM  AX 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 

h 

CALL 

CALL 

CALL 

CALL 

CALL 


=  AMT  Ml 
=  AMAXI 
PLOTS  (6 
RED 

PLOT  (l 
LAXTS  (9 
LAXIS  (7 
LAXTS  (0 

LOG LOG 
LOG LOG 
SYMBOL 
SYMBOL 
SYMBOL 


(PPy IN, PTHEOR (l ) , PTHEOR (NPOINT) ) 
(PPMAX, PTMAX) 


.,-6. 5,-3) 

.  , 0. ,  'PS  1**2 /HZ ’,9,6., 90. , PM AX , PM  IN , PDM I N , PDELTA ) 

. , 0. , 'PS  I **2/HZ ' ,-9,6. ,90. ,PMAX,PM IN, PDM IN, PDELTA) 

. ,0. , 'NORMALIZED  FREQUENCY' ,-20,7. ,0. ,W(NPOrNT) , 

W(l ) ,WDM IN , WDELTA ) 

(W, PTHEOR, NPOINT,  l , WDM  IN , WDELTA , PDM IN, PDELTA , 0,  '  . ’ ) 
(W, PPR IME , 102,1, WDM  IN, WDELTA, PDM IN, PDELTA,-!  ,  '* ' ) 

(1 . , 6. 2 , 0 , 'CHANNEL  10'  ,0. ,13) 

(3, , 6. 2,0, LABEL, 0. ,6) 

(3.75, 6,2, 0, 'SECONDS  AFTER  IGNITION ', , 22 } 
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0159 


PRINT  62 S 


01f»O 


0161 
0152 
0153 
D154 
0  1  6  5 

0  1 6  5 
0167 
0158 
0169 
0170 
0  171 
0  1  7  ? 
■0173 
017/1 
0175 
0176 

0177 

0178 


S'* 

V_. 

c 

A 

c 


950 


C 


1000 

c 


GOTO  1000 

CONVERT  W  INTO  REGULAR  FREQUENCY  AND  PLOT  ABOVE  CURVES  ON  LINEAR 
AXES. 


DO  950  1=1 , 102 
W(I)  =  W(I)*W0 
WINC  =  (W( 102 )  -  W(l))/7. 

CALL  PLOT  (0 . , -6 . , -3 ) 

PM  IN  =  AM1N1  (PPMIN ,  PTMIN ) 

PMAX  IS  THE  SAME  AS  ABOVE,  SO: 

PINC  =  (PMAX  -  PMIN)/6. 

CALL  AXIS1  (0. , 0. , 'PSI **2/HZ ' ,9,6. ,90. ,FMIN,PINC) 

CALL  AXIS1  (7 . , 0, , ' PSI **2 /HZ  1 ,-9,6. ,90. , PM  IN , PINC) 

CALL  AXIS1  (0. ,0. , 'FREQUENCY' ,-9,7. ,0. ,W(l ), WINC) 

CALL  LINE  (W,W(1 ), WINC, PTHEOR , PM  IN , PINC , 102,1,0, ) 
CALL  LINE  (W,W(i ), WINC, PPRIME,PMIN, PINC, 102, 1,-1, '*') 
CALL  SYMBOL  (1 . , 6. 2, 0, 'CHANNEL  10',0.,10) 

CALL  SYMBOL  (3 . , 6. 2 , 0 , LABEL , 0. , 5 ) 

CALL  SYMBOL  (3 . 75 , 6. 2 , 0 , ' SECONDS  AFTER  IGNITION ’,  0. , 22 ) 
PRINT  625 
CALL  CLOSE  (6 ) 

STOP 

END 
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C 


0001 
0002 
0003 
0004 
3005 
0006 
3  007 

0038  100 

0009 
00]  0 
001  1 


THIS  FUNCTION  RETURNS  THE  BEST  P-NOT  VALUE 
OMEGA-NOT.  CRITERION  IS  LEAST  SSE. 


FOR  A  GIVr 


FUNCTION  PNOT(WNOT) 

COMMON  PI,  W(300) ,  PPRIME  (2  56) 

DENSUM  =  9. 

TOPS  UM  =  0. 

DO  100  1=1,102 

WS'JM  =  W ( I  )  /WNOT  +  WNOT/W(  I  ) 

DENSUM  =  DENSUM  +  1./{WSUM**4) 

TOPSUM  =  TOPSUM  +  PPRIME (I )/ (WSUM* WSUM ) 
PNOT  =  PI  *  WNOT  *  TOPSUM/ (4.  *  DENSUM) 
RETURN 
END 
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000] 
0002 
000  3 
000/} 
0  00  5 


C 

C 

n 

c 

c 


THIS  FUNCTION  RETURNS  A  P (I }  VALUE  FOR  A  GIVEN  W(I) 
W(’  )  IN  OMEGA  FORM. 

FUNCTION  PCURVE (WNOT, P0,WI ) 

COMMON  PI 

WSUM  =  W I /WNOT  f  WNOT/WI 

PCURVE  =  4,  *  PH/ (PI  *  WNOT  *  WSUM  *  WSUM) 

RETURN 

F.ND 
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C 

C 

(KM  1 
'KM? 

9  Of  0 
0.0  0/1 
0  0  0  o 


0,i. in 


SUBROUTINE  FUNCT (VAR, WORK) 
DIMENSION  WORK (102) 

COMMON  P I , W ( 3 0 0 ) 

PO  =  PNOT  (VAR  ) 

DO  100  1=1,102 

W  )HK  ( l  )  =  PCURVE  (VAR,  P0,W(  I  )  ) 
RF.TURN 
ENO 
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0001 


SUBROUTINE,  FFTFP  (XREAL,  XIMAG  ,N,  M,  IF) 


I F=0  FORWARD  TRANSFORM 
I F=1  INVERSE  TRANSFORM 

M=0  XREAL  AND  XIMAG  RETURNED  AS  REAL  AND  IMAG.  FOR  FORWARD  X FOR MS 
M=1  "  "  "  "  "  MAGNITUDE  AND  PHASE  " 

(PHASE  IN  DEGREES) 

M=2  XREAL  RETURNED  AS  'PSD’;  XIMAG  =0. 

HERE  ' PS D 1  MEANS  SUM  OF  N  VALUES  OF  XREAL  =  MEAN  SQUARE  OF  INPUT 

FOR  INVERSE  TRANSFORM  M  DEFINITIONS  APPLY  TO  INPUT  DATA 
XREAL  AND  XIMAG  RETURNED  AS  REAL  AND  IMAGINARY 

FOR  FORWARD  TRANSFORMS  XREAL  AND  XIMAG  INPUT  AS  REAL  AND  IMAGINARY 


0002 
0003 
0004 
000  5 


DIMENSION  XREAL  (1  )  ,  XIMAG  (1  ) 
PI-3. 1415° 

DTOR=P 1/130. 

IF ( I F. EQ. 0 )GO  TO  6 


MUST  PREPARE  FOR  INVERSE  TRANSFORM 


0007 

0O09 


IF(M.EQ.0)GO  TO  2 
IF ( M. EQ. 2 )G0  TO  4 


0011 
0012 
001  3 
0014 
0015 
0015 


INPUT  IS  MAGNITUDE  AND  PHASE 

DO  1  1=1 ,N 
FMAG=XREAL  (I  }/N 

XREAL  (I ) =FMAG *COS (XIMAG (I )*DTOR) 
XIMAG  (I  )  =-FMAG*S  IN  (XIMAG  (I  )  *DTOR ) 
CONTINUE 
GO  TO  5 


INPUT  IS  REAL  AND  IMAGINARY 


0017 
0013 
0019 
0020 
00  21 


DO  3  1=1, N 

XREAL  (I  )=XREAL  (I  )/N 

XIMAG  (T  )=-XIMAG  (I  )/N 

CONTINUE 

GO  TO  5 


C 

0022  4 

0023 
00  24 
0025  5 


INPUT  IS  ’PSD' 

FACT=F  LOAT  (N  )  *N 
DO  5  1=1 ,N 

XREAL  (I  )  =XREAL  (I)/FACT 
CONTINUE 


00  2  5  5 

C 

0027 


CALL  FFTB  (XRF.AG ,  XIMAG  ,  N ) 
IF  ( I  F.  EQ.  1  ) R  -TURN 
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C 

C 


0O29 

0031 


00  3  3 
00  34 
0035 
00  39 
0037 
0033 


0039 
004  (’ 
9041 
1)042 
;K)43 
904  4 
904  5 


C 

C 

C 

3 


TRANSFORM  WAS  FORWARD 

IF ( M. EQ. 0) RETURN 
IF ( M. EQ. 2 )G0  TO  3 

DESIRE  OUTPUT  IN  MAGNITUDE  AND  PHASE 
DO  7  1=1 , N 

XMAG  =SQRT  (XREAL  (I  )  *XREAL  (I  )  FXIMAG  (I  )  *XIM  AG  (I  ) ) 

XIMAG  (I  )=ATAN2  (XIMAG  (I )  , XREAL  (I  ) )  /DTOR 

XREAL  (I )=XMAC 

CONTINUE 

RETURN 

DESIRE  'PSD' 

FACT=FLOAT  (N  )  *N 
DO  9  1=1, N 

XREAL  (I  )=  (XREAL  (I  )  *XREAL  (I  )+XIMAG  (I  )  *XIMAG  (I  )  )/FACT 
XIMAG  (I ) =9. 

CONTINUE 

RETURN 

END 
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PAG ft  001 


0001 
0002 
000  2 
0004 
000  5 
0005 
0007 

O00S  102 
0O09 
0010 
0011 
‘  0012 
001  3 
0014 
0015 
0015 
0017 
0018 
0019 
0020 

0021  lOi 
0822 
002? 

0025 

0025 

0027  100 

0028 

0029 

0  030 

0032 

O033 

00  34 

0035 

0035 

0  037 

0038  103 

0.0  39 
O04O 


SUBROUTINE  FFTB  (XRL‘AL,XIMAG,N) 
DIMENSION  XREAL  (N  )  ,  XIMAG  (N  ) 
NU=LOG2(N) 

N2=N/2 

NU1=N'J-1 

K-0 

DO  100  L=1,NU 
DO  101  1=1 ,N2 
P=IB  ITR  (K/2**NU1  ,NU) 

ARG =6 .28  3 185  *P/FLOAT (N ) 

C=COS  (ARG ) 

S=SIN  (ARG) 

Kl=K+l 
KIN  2=K l +N2 

TREAL=XREAL  (K  1N2  )  *C+XIMAG  (KIN 2  )*S 
TIMAG=XIMAG  (K  1N2)*C-XR£AL  (K1N2)*S 
XREAL  (!<  IN  2  )  =XREAL  (K 1  )  -TREAL 
XIMAG  (K1N2)=XIMAG  (K1  )-TIMAG 
XREAL  (K 1  )=XREAL  (K  1  )+TREAL 
XIMAG  (XI  )=XIMAG  (K1)+TIMAG 
K=K  +] 

K=K+N2 

IP (X . LT. N )GO  TO  10? 

K=3 

N'J  1  =NU1  - 1 

N2=N2/2 

DO  103  K=1,N 

I -IB  ITR  (K-l  f  N'J  )  +  l 

IF ( I . LE. K)GO  TO  103 

TREAL=XREAL  (K  ) 

TIM  AG  =X  I  MAG  (X) 

XREAL  (K)=XREAL  (I  ) 

XIMAG  (K)=XIMAG(t) 

XREAL  (I  )  =TREAL 

XIM  AG  ( I  )  =T IMAG 

CONTINUE 

RETURN 

END 
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OO  11 
0012  2 
0013  1000 

0014 
001  5 


Fri  09-Nov-79  00:00:oo 


FUNCTION  LOG 2  (N  ) 

N1=n 

■7=1 

LOG 2=0 

IFn.EQ.Mi  )return 
fF(  7.GT.N1  )00  TO  7 
■7=7*2 

LOG2=LOC2+l 

GO  TO  l 
TYPE  100O, nj 

format  (IX. 15,'  IS  NOT  A  POWER  0p  ,,, 
END 
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V  .  0  'I  A 

SUBROUTINE  HANN  (R  IN ,  N  ) 

0002 

DIMENSION  RIN (1 ), TEMP (256) 

0003 

A  =M  - 1 

••)”,/ 4 

00  1  I  ~2  , M 

,5  Oil  5 

TEMP  (I  )=  (RIN  ( I -1  )  l-RIM  ( f  +1  ))/4+RIN(I  )/2 

?!  0  0  A  I 

CONTINUE 

fitlU? 

RIN  (1  )  =  (R  I N  (l  )  HR  IN  (2  )  )/2 

0008 

RTN  (N)=  (RIN  (N)  *-RIN  (M)  ) /2 

0009 

DO  2  I  =2  ,  M 

0010 

RIN  (I  )=TEMP  (I  ) 

0011  2 

CONTINUE 

0012 

RETURN 

0013 

end 
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no  01 
0002 
0033 
0  03  4 
3035 
0035 

no  07 
0033 
0009 
3010 
001  1 
1)01.2 
0‘".  13 


SUBROUTINE  HAN 2  (R IN , N ) 

DIMENSION  RIN  (l  )  , TEMP  (2  55) 

M=N  -1 

00  1  1=7,  M 

TEMP  (I  )=  (RIN  (1-1  )+RIN  (T  +1  )  1/3.+RIN  (I  )/3, 
CONTINUE 

RIN  (1  )=  (RIN  (1  )*RIN  (2  )+RIN  (3  )  }/7 
RIN  (N)=  (RIN  (N  )  +  RIN  (M)  +-RIN  (M-l  )  )/3 
DO  2  1=2,  M 
RIN  (I  )=TEMP(I  ) 

CONTINUE 

RETURN 

END 
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C 
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f* 

C 
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r* 

r 

n 

r 

r 

r 

r 

r 

i' 
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0 
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r 
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C 


SUBROUTINE  SEARCH (THETA , SSETHE , AMPTHE , IAMP , ALOW , HIGH , ACCTHE 
I ACCSE , WORK , STAN , NPTS / I PR I NT, I  PLOT, X3 12 E , YS IZE , VARNAM , NCHAR 
VARPLT, SSEPLT, NPLTPT ) 

PURPOSE 

DOES  A  "STAB"  SEARCH  FOR  AN  UNKNOWN  PARAMETER. 


ARGUMENTS 


THETA 

VALUE  OF  PARAMETER  RETURNED  AFTER  THE  SEARCH. 
SSETHE 

SSE  ASSOCIATED  WITH  THETA,  RETURNED  BY  PROGRAM. 
AMPTHE 

AMPLITUDE  RATIO  FOUND  AND  RETURNED.  IF  NO  RATIO 
CALCULATIONS  ARE  DESIRED,  SET  IAMP  =  0 

IAMP 

SEE  ABOVE. 

ALOW 

USER-SUPPLIED.  LOW  END  OF  RANGE  IN  WHICH  THETA  TS 
EXPECTED  TO  APPEAR. 


HIGH 

USER-SUPPLIED.  HIGH  END  OF  RANGE  FOR  THETA. 

TF  THE  MINIMUM  SSE  IS  FOUND  AT  EITHER  ALOW  OR  HIGH, 

THE  PROGRAM  WILL  EXPAND  THE  RANGE  UNTIL  A  NON-BOUNDARY 
MINIMUM  SSE  IS  FOUND. 


ACCTHE 

USER-SUPPLIED.  THE  SEARCH  ENDS  WHEN  EITHER  THE 
PARAMETER  OR  ITS  SSE  HAS  BEEN  FOUND  TO  SPECIFIED 
PRECISIONS. 

WHEN  THE  DIFFERENCE  BETWEEN  SUCCESSIVE  ESTIMATES  OF  THE 
PARAMETER  TS  LESS  THAN  OR  EQUAL  TO  ACCTHE,  THE  SEARCH 
IS  STOPPED. 

RETTING  ACCTHE  =  3.  RESULTS  IN  THE  COMPUTER  SEARCHING  TO 
1’HE  LIMITS  OF  ITS  OWN  ACCURACY. 

I ACCSE 

WHEN  SUCCESSIVE  SSE '5  AGREE  TO  I ACCSE  DIGITS  THE  SEARCH 
I S  STOPPED.  MAXIMUM  =  7 

WORK 

WORK  VECTOR.  LENGTH  =  NPTS. 


P 


43 


FORTRAN 


VO 2.; M 


Fri  09-No v-79  00:00:00 


PAGE  00  0 


STAN 

USER -SUPPLIED  DATA  VECTOR  CONTAINING  THE  POINTS  THF 
PROGRAM  TRIES  TO  MATCH. 

NPTS 

NUMBER  OF  POINTS  IN  STAN. 

I PR I NT 

UNLESS  SET  =  0,  USER  GETS  PRINTOUT  OF  INTERMEDIATE 
RESULTS  OF  THE  SEARCH. 

I  PLOT 

PLOTTING  CONTROL  CHARACTER 
=0  NO  PLOT 

=1  PLOT  WITH  REGULAR  AXES 

=2  SEMI  LOG  PLOT  (SSE'S  ON  LOG  AXIS) 

=1  30TH  PLOTS 


XSIZE 

SIZE  OF  X-AXIS  OF  THE  OPTIONAL  PLOT. 

YSIZE 

SIZE  OF  Y-AXIS  OF  THE  OPTIONAL  PLOT. 

'YARN  AM 

HOLLERITH  STRING  NAME  OF  PARAMETER. 

MCHAR 

NUMBER  OF  CHARACTERS  IN  VARNAM. 

VARPLT 

ARRAY  USED  FOR  PLOTTING.  LENGTH  =  NPTPLT 

UPON  RETURN  CONTAINS  ALL  ESTIMATES  OF  THE  VARIABLE, 

ORDERED  SMALL  TO  LARGE. 

G5EPLT 

ARRAY  USED  FOR  PLOTTING.  LENGTH  =  NPTPLT 

UPON  RETURN  CONTAINS  SSE  FOR  EACH  VARIABLE  ESTIMATE 

NPTPLT 

SIZE  OF  PLOTTING  ARRAYS. 

IF  THE  INITIAL  RANGE  (HIGH  -  ALOW)  HOLDS  THE  BEST  VALUE 
OF  THE  UNKNOWN  VARIABLE,  THEN  THE  FOLLOWING  FORMULA  SHOULD 
PROVIDE  AN  UPPER  BOUND  ON  NPTPLT: 

NPTPLT  <=  3  +  2 * (LOG (RANGE/ ACCTHE ) ) /LOG  (2 ) 

THE  USER  MUST  WRITE  A  SUBROUTINE  FUNCT (VAR, WORK)  THAT  GENERATES 
A  VECTOR  "WORK"  OF  ESTIMATED  DATA  POINTS,  GIVEN  THAT  THE  UNKNOWN 
PARAMETER  EQUALS  "VAR"  . 


-SUBROUTINE  SEARCH  (THETA ,  SSETHE ,  AMP7HE ,  I A.MP,  ALOW,  HIGH ,  ACCTHF. , 

I ACCSE , WORK , STAN , NPTS , I PR  I NT , I PLOT , XS I Z  E , YS I Z  E , VARNAM , NCHAR , 
i-  VARPLT,  SSEPLT,NPLTPT) 
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PACK  uoi 


11  ,1  .1  .> 

,n 

1 1 '  i » *  •* 


dimension  vmm'O  ,bsf.(A)  .amima)  ,work(npts)  ,stan  (npts)  , 

VARPLT  (NPLTPT)  ,  SSEPLT  (NPLTPT) 

LOGICAL*  1  ACCHE  1  (14)  ,  ACCSE?  (14)  ,  VARNAM  (N01IAR )  ,  ALPH  (  ’  7  ) 

OAT  A  ALPH/  ’  *  ,  *  A  ’  ,  •  ,  *C  ’  ,  *D  '  ,  '  E  '  ,  '  F  '  ,  'G  •  ,  '  11  •  ,  '  1  '  ,  '  1  •  ,  ’K  '  , 

'L  1  ,  'M'  ,  'N  '  ,  'O'  ,  'P  '  ,  '0'  ,  'R  '  ,  *S  '  ,  *T'  ,  M)  '  ,  'V  '  ,  1  W  ,  ’X  '  ,  'V  1  ,  '7.  '/ 

1K(1A*TBF.  .GT.  '> )  STOP  ' G H L  ACCURACY  CONSTANT  MUST  PR  <o  /' 


O  l  FOLD  0.  !  USED  IN  PRECISION  CHECK  SECTION 
1  ROUND  a  1  !  MARKS  THE  ROUND  JUST  COMPLETED 

NPOINT  a  0  !  MARKS  THE  MUMPER  OK  POINT1'  ENTERED  INTO 
THE  PLOTTING  ARRAYS 

1  STOP  »  0  !  IK  THE  PLOTTING  ARRAYS  ARE  KILLED,  SSRSIH 

SUP  ROUTINE  SETS  THIS  EQUAL  TO  1,  WHICH  1M 
TURN  CAUSES  THE  SEARCH  TO  STOP. 

NOTE  l  !  INDEX  OK  THE  "ALPH"  ARRAY.  USED  WHEN  THE 
SEARCH  RANGE  IS  EXTENDED. 

ADD  (HIGH  -  A  LOW ) *  .  7  A 


Ov!  1  1 

Ui'-l  4 


s'  1  A 
UP  l  A 


,LI  1  ;1 


Pi!,”' 

Pit  M 
■ip.’  >  l 


>ET  U"  INITIAL  VAR  ( 1  )  ,  VAR  (I  ) ,  AND  VAR  (A  )  WITH  GRE'S 
VAR  (A)  -  HIGH 

CALL  SSKGUP  (VAR  (A  )  ,  SHE  (A  1  ,  AMP  (A  )  ,  1  AMP, WORK  ,  STAN  ,  NPTS  ,  1  PL  U', 
1  NPO l NT, VARPLT, SSEPLT, NPLTPT, 1  STOP ) 

VAR  (  1  1  *  ALOW  t  .  A*  (11  rC.ll-ALOW) 

CALL  uSESUP  (VAR  (  !  )  ,SSE  (  * )  ,AMP  (.4  )  ,  1  AMP , WORK ,  STAN  ,  NPTS ,  1  PLOT, 
i  NPOl  NT,  VARPLT,  SSEPLT,  NPLTPT,  1  ST D I' ) 


OP*  1  ir.l  V  MM  l  )  --  PAR  (  l  )  (VAR  (A  )  -  VAR  (  1  )  ) 


1  F  (VAR  (  1 )  .NE.  VAR  ( l  )  )  GOTO  i:>A  (PRECISION  CHECK  REQUIRED  PECAU*  '• 

RANGE  EXTENSION  HF NDS  CONTR'l. 

HACK  HERE. 

1  START  *  '> 

GOTO  4  14 

CALI.  SSf.SUH  (VAR  (1  )  ,SSE(1  )  ,AMP(1  )  ,  l  AMP ,  WORK  ,  STAN  ,  NPTS  ,  1P1.0T, 

•  NPOINT,  VARPLT,  SSEPLT ,  NPLTPT ,  !  s-pOP  ) 


c  '-ET  UP  VAR  CM  WITH  SHE  (2  ) 

r 

1  VAR  CM  *  (VAR  (  1  )  -  VAR  (l  ))*.*>  1  VAR(l) 

CALL  SSEHU’.l  (VAR  C>  )  ,  SSE  ('? )  ,AMP  (?.  )  ,  TAMP , WORK ,  STAN ,  NPTS  ,  I  PLOT, 
i  NPOl NT, VARPLT, SSEPLT, NPLTPT , l STOP ) 

< ' 

o  IK  EITHER  OK  SHE  (l  )  OR  SHE  (2 )  IS  THE  MINIMUM,  THEN  THERE  IS  NO 

r  NEED  TO  FIND  SHF,  (4  )  .  IF  SHE  ( 1  )  IS  MIN,  EXTEND  THE  RANGE  DOWNWARD,. 

r  IK  SSL  (  M  IS  MIN,  GO  ON  TO  PRECISION  CHECKS  BEFORE  DOIN'’,  ANOTHER 

c  ROUND.  OTHERWISE,  COMPUTE  VALUES  FOR  VAR(1). 

<  * 

IF((GSM1)  .GT.  HSE(M)  .OR.  (SHF.  (l)  .GT.  SHE  (  M  )  .OR. 
i  (SSE  ( 1  )  .GT.  SHE  (A)))  GOTO  MU 

»K(1 PRINT  .EQ.  -M  GOTO  I/O 
"HINT  ImO,  1  ROUND,  ALPH  (Noti-M  ,  VAR,  SSE 
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093') 
P'.Pl 
0032 
0  0  3  4 


0  0  $7 
0033 
0  0  3  0 
'  )  9  4  \'l 


FORMAT  ('  ROUND',  [3, Al  /  '  VAR  I  ABLE  :  1  5G  1  5 . 5  /  '  SSE:'5X,  *'  J  1  -  i  ! 
NOTE  =  MOVE  +  l 
7 "  (NOTE  ,£Q.  23)  NOTE  =  27 
•*  L  CLOSE  (3  ) 

IF  (I  STOP  .NE.  0)  GOTO  #550 

RESET  ARRAYS  AND  DO  ANOTHER  ROUND: 

CALL  RESET  (VAR, 1,3, 4, 5) 

CALL  RESET (SSE, 1,3, 4, 5) 

CALL  RESET (AMP, l , 3, 4, 3) 


00  4  3 
0.04  4 


I  F  (  (SSE  (2  }  .GT.  SSE  (3))  .OR.  (SSE  (2  )  .GT.  SSE  (3  )  )  )  GOTO  Is  4 
HERE  IT  IS  KNOWN  THAT  VAR  (2 )  IS  THE  BEST  VALUE  SO  FAR 
r START  =  \ 

GOTO  420 


00  4  3  30v) 
0  049 


0O4  7 
004  9 
00  51 
0  052 
0  0  53 
00  55 
0  055 
0058 
0059 
0  0  59 
00  51 
0052 
(10  54 
■0  0  55 


0  057 


SET  UP  VAR  (4)  WITH  SSE  (4  ) 

VAR  (4)  =  (VAR  (5)  -  VAR  (3))  *.5  +  VAR  (3) 

CALL  SSESUB (VAR  (4 ) ,SSE  (4 ) ,AMP  (4 ) , I AMP, WORK, STAN, NPTS,  I  PLOT, 

-  NPOINT, VARPLT, SSEPLT , NPLTPT , 1  STOP ) 

AS  ABOVE,  SEE  IF  THE  ENDPOINT  OF  THE  SEARCH  INTERVAL  (HERE,  V\«  ('• 
YIELDS  A  BOUNDARY  MIN  SSE.  IF  SO,  RESET  THE  VARIABLES  SO  THAT 
THE  RANGE  IS  EXTENDED  UPWARDS.  OTHERWISE,  FIND  THE  STARTING  PC’.  N 
FOR  THE  NEXT  ROUND,  AND  GO  TO  PRECISION  CHECK  SECTION. 

IF  (  (SSE  (5)  .GT.  SSE  (3))  .OR.  (SSE  (5)  .GT.  SSE  (4))  )  GOTO  409 

IF (I  PRINT  . EO.  0)  GOTO  350 

PRINT  150,  TROUNO,ALPH  (NOTE) , VAR, SSE 

NOTE  =  NOTE  t-  1 

IF (NOTE  .EO.  23)  NOTE  =  27 

CALL  CLOSE  (5) 

IF (I  STOP  .NE.  0)  GOTO  550 
CALL  RESET (VAR, 3, 1,2, 3) 

CALL  RESET (SSE, 3, 1,2,  3) 

CALL  RESET  (AMP, 3,1,2, 3) 

VAR  (5)  =  VAR  (3)  +  (VAR  (3)  -  VAR  (1  )  ) 

IF  (VAR  (5)  .NE.  VAR  (3))  GOTO  375 
I START  =  2 
GOTO  434 

CALL  SSESUB  (VAR  (5 ) ,SS5  (5) ,AMP(5) , I AMP , WORK , STAN , NPTS  ,  T  PLOT  . 
i-  NPOINT,  VARPLT,  SSEPLT,  NPLTPT,  I  STOP) 

GOTO  300 


0058  400 

0  0  59 

9071  420 

0073 
0R7  4 

C 


l  START  2 

IF  (SSE  (4  )  .  LE .  SSE  (3))  ISTART  =  3 

IF ([PRINT  ,EQ.  0)  GOTO  430 

PRINT  150,  I  ROUND, ALPH  (NOTE) , VAR, SSE 

CALL  CLOSE  (5) 

PR  EC  1  ;  s  ON  CHECKS 


F I RST : 


IF  THE  DIFFERENCE  BETWEEN  GUCCF.G  1 VF.  ESTIMATES  1  !' 
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r 

C 

C 

#*» 

C 

c 

c 

0075  430 

,107  5 
i)079 

(5079  V'l 

0  c>  3  n 

0031  432 

0033  434 

003  <1  4  35 


0  U  3  5 


UNKNOWN  VARIABLE  IS  LESS  THAN  OR  EQUAL  TO  " ACCTHE" 
(USER-SUPPLIED  CONSTANT)  THEN  EXIT. 

THE  DIFFERENCE  MAY  NEVER  BE  LESS  THAN  ACCTHE  BECAUSE  OF 
ROUNDING,  IN  THIS  CASE  DIFF  REMAINS  CONSTANT  OVER  ROUNDS. 
DIFOLD  IS  A  CHECK  ON  THIS. 

DIFF  =  VAR  (2  )  -  VAR  (I ) 

1 F (D IFF  .GT.  ACCTHE)  GOTO  432 
PRINT  431, ACCTHE 

FORMAT  (’0  *****  SEARCH  STOPPED  DIFFERENCE  BETWEEN  ESTIMATES  <  = 

+  G15.5) 

GOTO  700 

IF (D IFF  .NE.  DIFOLD)  GOTO  440 
PRINT  435, VARNAM 

FORMAT ( '0  *****SEARCH  STOPPED*****  '  / 

+  '  DOUBLE  PRECISION  REQUIRED  FOR  FURTHER  REFINEMENTS  OF', 

4-  *  ESTIMATES  OF  ',30Al) 

GOTO  700 


C 

C 

r 


SECOND:  IF  THE  SSE'S  AGREE  TO  A  USER-SUPPLIED  NUMBER  OF 

DIGITS  (IACCSE)  THEN  EXIT. 


00  96 

r 

44  3 

DIFOLD  =  DIFF 

0  037 

ENCODE  (14,450, ACCSE 1 ( 1 ) ) SSE (I  START ) 

’.5038 

4  50 

FORMAT (1  PE  14.7) 

■HJ39 

DO  5  50  1=  (ISTART+l )  , (ISTART+2) 

09  9(5 

ENCODE  (14,450, ACCSE2  (i  )  )SSE  (I  ) 

0  09  3 

DO  500  1=12,14 

0  092 

IF  ( ACCSE  1  (I  )  .NE.  ACCSE2  (J  )  )  GOTO 

503 

U094 

50  0 

CONTINUE 

0.-95 

DO  550  J=2 , ( IACCSE  +  2) 

00 

I F  (  ACCSE  1(1)  .NE.  ACCSE2H))  GOTO 

500 

0  09  9 

5 

CONTINUE 

9  499 

PRINT  575, [ACCSE 

o ;  4 

57  5 

FORMAT  (’l)  *****  SEARCH  STOPPED  SUM 

^  L  2  '  DIGITS.') 

OF  SQUARED 

0  1  01 

GOTO  700 

•4ET  UP  ARRAYS  FOR  THE  NEXT  ROUND 

4  i  02 

■V4'4 

[ F ( I STOP  .NE.  C)  GOTO  550 

9  1  0  4 

CALL  RESET  (VAR,  I START,  1,3,5) 

0  1  9  5 

CALL  RESET (SSE, I  START , 1 , 3 , 5 ) 

9  I  0 .9 

CALL  RESET (AMP, I START , 1,3.5) 

'4  I  47 

[ROUND  =  [ROUND  +  1 

9  i  0  8 

NOTE  =  1 

01  09 

GOTO  1 50 

9  !  !  0 

4  50 

PRINT*,'*****  SEARCH  STOPPED  PLOTTING  ARRAYS 

0  ’.  !  1 

7  0!) 

THETA  =  VAR ( f  START +  1  ) 

0  1  1  2 

GOETHE  =  SSE  ( [ START+1 ) 

0. 1  i  ; 

C 

AMRTHE  =  AMP (i START+1  ) 

o  i  i 

L  F { I  PLOT  .EQ.  fi)  RETURN 

MATCH  TO' 


CO  ****** 
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0  3  15 
0117 
0118 
,  0119 
0120 
0122 
0123 
‘  0124 
0125 
0125 
0127 
0128 

0129 

01  31 

0132 


0133 

0135 

0136 

1  37 
0133 


C 

C 

C 

C 

830 


850 


900 

C 


C 


C 


C 

C 

950 

C 


ORDER  THE  VARIABLE  ESTIMATES  SO  THAT  A  LINE  CAN  BE  DRAWN  BETWEEN 
DATA  POINTS  ON  THE  GRAPH,  SHOWING  BEHAVIOR  OF  SSE'S. 

FORMAT  ( ’ 1 ' ) 

DO  903  1=1 , (NPOTNT  -  1) 

MIN  =  I 

DO  850  J=I , NPOINT 

IF (VARPLT (3 )  .LT.  VARPLT (MIN ) )  MIN  =  J 
CONTINUE 
SI  =  VARPLT  (I) 

VARPLT  (I)  =  VARPLT (MIN) 

VARPLT  (MIN)  =  SI 
SI  =  SSEPLT  (I  ) 

SSEPLT(T)  =  SSEPLT(MIN) 

SSEPLT (MIN )  =  SI 

IF ( I  PLOT  .EQ.  2)  GOTO  950 

PRINT  800 

REGULAR  PLOT  OF  SSE  VS  VARIABLE  ESTIMATES. 

CALL  DOPLOT  (XS IZ E , YS IZ E , VARPLT, SSEPLT, NPOINT, VARNAM , NCHAR , ’SSE ' , 

3,0,  '  .',1,1) 

IF ( I  PLOT  .EQ.  1)  RETURN 

PRINT  800 

SEMILOG  PLOT  OF  SSE  VS  VARIABLE  ESTIMATES. 

CALL  PLTLGY  (XS IZ E , YSIZ E , VARPLT, SSEPLT , NPOINT, VARNAM , NCHAR , 'SSE '  , 

3,0, '.\n 

RETURN 

END 


t 

3 


4 

i 


i 
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0001 

0002 

0033 

000  4 
000  5 

0  007 
000  3 
0009 
0010 
0011 
0012 
0  013 


09  M 
9015 
0010 
0017 


90  13 
0  0  20 
0  02  1 
0022 
0023 


O0  2  5 
90  2  5 


C 

C  SUBROUTINE  SSESUB 

C 

C  PURPOSE 

C  CALLED  BY  SEARCH  TO  FIND  AMPLITUDE  RATIO  (AMP)  AND  SSE  FOR  A 

C  GIVEN  VALUE  OF  THE  UNKNOWN  PARAMETER  (VAR ) 

C 

SUBROUTINE  SSESUB (VAR , SSE , AMP , I AMP , WORK , STAN , NPTS , I  PLOT , 

+  NPOINT, VARPLT, SSEPLT, NPLTPT, ISTOP) 

DIMENSION  WORK  (NPTS)  ,  STAN  (NPTS)  ,  VARPLT  (NPLTPT)  ,  SSEPLT  (NPLTPT) 

C 

CALL  FUNCT  (VAR, WORK) 

C 

A  =  l. 

I F ( I  AMP  .EQ.  0.)  GOTO  200 
C  FIND  AMPLITUDE  RATIO 

SUMSQ  =  0. 

SUMX  =  0. 

DO  100  1=1, NPTS 

SUMX  =  SUMX  +  STAN  (I  )*WORK  (I  ) 

100  SUMSQ  =  SUMSQ  +  WORK (I ) *WORK  (I ) 

AMP  =  3UMX/SUMSQ 
A  =  AMP 
C 

C  FIND  SSE 

20w  SSE  =  0. 

DO  300  1=1, NPTS 

ERROR  =  STAN  (I)  -  A*WORK(I) 

300  SSE  =  SSE  +  ERROR*  *2 

C 

C  FILL  PLOTTING  ARRAYS  IF  NECESSARY 

C 

r  F ( T PLOT  .EQ.  o)  RETURN 
NPOINT  =  NPOINT  +  I 
VARPLT  (NPOINT)  =  VAR 
SSEPLT  (NPOINT)  =  SSE 

■190  IF  (NPOINT  .GE.  (NPLTPT  -  1)  )  ISTOP  =  1 

C  IF  ISTOP  =  1,  ARRAYS  ARE  FILLED  AND  SEARCH  STOPS  UPON  RETURN. 

C 

RETURN 

END 
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r* 

c 

c 

c 

c 

0001 
0002 
000  3 
0304 
00  05 
000  6 
0037 
0008 
0309 
0310 


RAG 


SUBROUTINE  RESET 
PURPOSE 

CALLED  BY  SEARCH  TO  REASSIGN  ELEMENTS  OF  ARRAYS 

SUBROUTINE  RESET (A, I , J, K, L) 

DIMENSION  A ( 5 ) 

X  =  .*.  { I ) 

Y  =  A  (  T  +1  ) 

Z  =  A  ( I  +2  ) 

A(J)  =  X 
A(K)  »  Y 
A  ( L  )  =  7. 

RETURN 

END 
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